TECHNICAL  REPORT  ARCCB-TR-95006 


AD 


ADAPTIVE  FINITE  ELEMENT  METHOD  IV: 
MESH  MOVEMENT 


J.M.  COYLE 
J.E.  FLAHERTY 


FEBRUARY  1995 


US  ARMY  ARMAMENT  RESEARCH, 
DEVELOPMENT  AND  ENGINEERING  CENTER 

CLOSE  COMBAT  ARMAMENTS  CENTER 
BENET  LABORATORIES 
WATERVLIET,  N.Y.  12189-4060 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


19950425  056 


DISCLAIMER 


The  findings  in  this  report  are  not  to  be  construed  as  an  official 
Department  of  the  Army  position  unless  so  designated  by  other  authorized 
documents . 

The  use  of  trade  name(s)  and/or  manufacturer(s)  does  not  constitute 
an  official  indorsement  or  approval. 


DESTRUCTION  NOTICE 

For  classified  documents,  follow  the  procedures  in  DoD  5200. 22-M, 
Industrial  Security  Manual,  Section  11-19  or  DoD  5200. 1-R,  Information 
Security  Program  Regulation,  Chapter  IX. 

For  unclassified,  limited  documents,  destroy  by  any  method  that  will 
prevent  disclosure  of  contents  or  reconstruction  of  the  document. 

^  For  unclassified,  unlimited  documents,  destroy  when  the  report  is 
no  longer  needed.  Do  not  return  it  to  the  originator. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


PuDiic  feoonmg  burden  for  this  coneaion  of  information  is  estimated  to  average  1  hour  per  response,  incJud»ng  the  time  for  reviewing  mstruaicns.  searcr^mg  e*  stmg  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  ana  reviewing  the  coUection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  otr'er  aspea  of  this 
colleaion  of  information,  including  suggestions  for  reducing  this  burden,  to  i/Vashmgton  Headquarters  Services,  Directorate  for  information  Operations  and  fteports.  1215  Jefferson 
Davis  Highway,  Suite  1204.  Arlington.  VA  22202*4302.  and  to  the  Office  of  Management  and  Budget.  Paperwork  Reduaion  Project  (0704-0 1S8),  Washington,  DC  20503. 


1.  AGENCY  USE  ONLY  (Leave  blank) 


2.  REPORT  DATE 
Fcbruaiy  1995 


4.  TITLE  AND  SUBTITLE  _ 

ADAPTIVE  FINTTE  ELEMENT  MEmOD  IV: 

MESH  MOVEMENT 


6.  AUTHOR(S) 


REPORT  TYPE  AND  DATES  COVERED 
Fmal 


5.  FUNDING  NUMBERS 


J.M.  Coyle  and  J.E.  Flaheity 


AMCMS:  6126.24.H19L1 
PRON:  1A17Z1CBNMBJ 


7.  PERFORMING  ORGANIZATION  NAME(S)  ANO  AOORESS(ES) 
U.S;  Anny  ARDEC 

Bendt  Laboratories.  AMSTA-AR-CCB-O 
Water^l,  NY  12189-4050 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

ARCCB-TR-95006 


9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  ADOet^lKSl  for 


U.S.  Army  ARDEC 

Qose  Combat  Armaments  Center 

Picatinny  Arsenal,  NJ  07806-5000 


11.  SUPPLEMENTARY  NOTES 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 


Approved  for  public  release;  distribution  unlimited 


NTIS  CRA&I 
DTIC  TAB 
Unannounced 
Justification 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


. 

Distribution  / 


Avail  and/or 
Special 


12b.  DISTRIBUTION  CODE 


13.  ABSTRACT  (Maximum  200  words) 

An  adaptive  finite  element  method  is  developed  to  solve  initial  boundary  value  proUems  for  vector  systems  of  parabolic  partial 
dififerential  equations  in  one  space  dimension  and  time.  The  differential  equations  are  discretized  in  ^ace  using  piecewise  linear 
finite  element  approximations.  Superconveigence  properties  and  quadratic  polynomials  are  used  to  derive  a  computationally 
inexpensive  approximation  to  the  spatial  component  of  the  error.  This  technique  is  coupled  with  time  integration  schemes  of 
successively  hi^er  orders  to  obtain  an  approximation  of  the  temporal  and  tot^  discretization  errors.  Tbe  stability  of  several 
mesh  equidistribution  schemes  for  time-dependent  partial  differential  equations  is  studied.  The  sdiemes  move  a  finite  difference 
or  finite  element  mesh  so  that  a  given  quantity  is  uniform  over  the  domain.  Mesh  moving  methods  that  are  based  on  solving 
a  system  of  ordinary  differential  equations  for  the  mesh  velocities  are  considered  and  some  of  these  methods  are  shown  to  be 
unstaUe  with  respect  to  an  equidistributing  mesh  when  the  partial  differential  q/stem  is  dissip.ative.  Simple  criteria  for 
determining  the  stability  of  a  particular  method  are  developed  and  the  construction  of  staUe  differential  systems  for  the  mesh 
velocities  is  demonstrated.  Several  examples  illustrating  stable  and  unstaUe  mesh  motions  are  present 


^^Fara^ltc  £i§^ntial  Equations,  Adaptive  Finite  Elements,  Equidistribution, 
Linear  Stability,  Mesh  Movement 


15.  NUMBER  OF  PAGES 


16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION  118.  SECURITY  CLASSIFICATION  119.  SECURITY  CLASSIFICATION  I  20.  LIMITATION  OF  ABSTRACT 

OF  REPORT  OF  THIS  PAGE  I  OF  ABSTRACT  I 


UNCLASSIFIED 


NSN  7540.01-280-5500 


UNCLASSIFIED 


UNCLASSIFIED 


Standard  Form  298  (Rev  2-89) 

p’o^cT'oea  bv  Std  Z39*’8 


TABLE  OF  CONTENTS 


INTRODUCTION  .  1 

STATIC  SPATIAL  EQUIDISTRIBUTION  .  4 

Example  1  .  5 

MESH  DYNAMICS:  LINEAR  STABILITY .  8 

Example  2 .  9 

Example  3 .  9 

MESH  DYNAMICS:  NEUTRAL  STABILITY .  14 

Example  4 .  15 

Example  5 .  16 

Example  6 .  17 

Example  7 .  17 

MESH  MOVEMENT  IMPLEMENTATION .  23 

SUMMARY . 26 

REFERENCES .  27 

List  of  Illustrations 

1.  Mesh  trajectories  for  Example  2 .  11 

2.  Graphs  of  Eq.  (27)  for  selected  instances  of  time .  12 

3.  Mesh  trajectories  for  Example  3 .  13 

4.  Mesh  trajectories  for  Example  4 .  19 

5.  Mesh  trajectories  for  Example  5 .  20 

6.  Mesh  trajectories  for  Example  6 .  21 

7.  Mesh  trajectories  for  Example  7 . 22 


1 


INTRODUCTION 


This  is  the  fourth  in  a  series  of  four  reports  whose  overall  purpose  is  to  describe  an 
adaptive  finite  element  method  (AFEM)  for  solving  systems  of  parabolic  partial  differential 
equations.  In  particular,  AFEM  attempts  to  find  a  numerical  solution  of  an  M-dimensional 
system  of  the  form 

,  a<x<b  ,  f>0,  (la) 


subject  to  the  initial  conditions 

u(a:,0)  =  u\x)  ,  a^x^b 


and  linear  separated  boundary  conditions 

A\t)u(b,t)  +  B\t)uJ,b,t)  =  g\t)  ’ 


(lb) 


(Ic) 


The  variables  x  and  t  represent  spatial  and  temporal  coordinates  and  denote  partial 
differentiation  when  they  are  used  as  subscripts;  u,  f  u°,  and  g  are  A/-vectors;  and  D,  A\  B,  A% 
and  R  are  MxM  matrices. 

The  problem  is  assumed  to  be  well-posed  and  parabolic;  thus,  e.g.,  D(x,t,u)  is  positive 
definite.  The  rows  of  B  and  R  are  restricted  to  be  either  entirely  zero  or  a  row  of  the  MxM 
identity  matrix.  When  the  z"*  row  of  or  is  identically  zero,  then  A^  or  Afi  cannot  be  zero, 
respectively,  and  the  boundary  condition  is  a  Dirichlet  (essential)  condition.  Otherwise,  the 
boundary  condition  is  a  Neumann  or  Robbins  (natural)  condition.  The  ultimate  goal  of  AFEM 
is  to  determine  an  approximate  solution  to  Eq.  (1)  to  within  a  user  prescribed  error  tolerance. 

The  adaptive  strategies  utilized  by  AFEM  are  (1)  error  estimation  coupled  with  (2)  local 
mesh  refinement  (cf.,  e.g.,  Adjerid  and  Flaherty  (ref  1),  Babuska  and  Dorr  (ref  2),  Babuska, 
Zienkiewicz,  Gago,  and  Oliveira  (ref  3),  Bank  and  Weiser  (ref  4),  Berger  and  Oliger  (ref  5), 
Bieterman  and  Babuska  (refs  6,7),  Moore  and  Flaherty  (ref  8),  Shephard,  Qingxiang,  and 
Baehmann  (ref  9),  Strouboulis  and  Oden  (ref  10),  Zienkiewicz  and  Zhu  (ref  11)),  and  (3)  mesh 
movement  (cf.,  e.g.,  Adjerid  and  Flaherty  (ref  1),  Arney  and  Flaherty  (ref  12),  Bell  and  Shubin 
(ref  13),  Davis  and  Flaherty  (ref  14),  Dorfi  and  Drury  (ref  15),  Dwyer  (ref  16),  Ewing,  Russell, 
and  Wheeler  (ref  17),  Hyman  (ref  18),  Kansa,  Morgan,  and  Morris  (ref  19),  Miller  and  Miller 
(refs  20,21),  Petzold  (ref  22),  Rai  and  Anderson  (ref  23),  Russell  and  Ren  (ref  24),  Saltzman  and 
Brackbill  (ref  25),  Smooke  and  Koszykowski  (ref  26),  Thompson  (ref  27),  Verwer,  Blom, 
Furzeland,  and  Zegeling  (ref  28),  and  White  (ref  29)). 

The  purpose  of  this  report  is  to  describe  the  mesh  moving  procedures  employed  by 
AFEM.  Detailed  summaries  of  how  AFEM  implements  its  other  adaptive  strategies  are  found  in 
separate  reports  entitled:  Adaptive  Finite  Element  Method  II:  Error  ^timation  (ref  30)  and 
Adaptive  Finite  Element  Method  III:  Mesh  Refinement  (ref  31).  Furthermore,  Adapative  Finite 
Element  Method  I:  Solution  Algorithm  and  Computational  Examples  (ref  32)  describes  how  all 
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the  adaptive  algorithms  are  implemented  in  unison  and  contains  results  demonstrating  the  utility 
of  AFEM  as  a  computational  tool. 

The  mesh  movement  performed  by  AFEM  is  based  on  equidistribution.  An 
equidistributed  mesh  is  a  partition  of  a  spatial  domain  into  elements  such  that  a  prescribed 
quantity  is  uniform  over  the  mesh.  More  specifically,  given  an  interval  (a,b)  and  a  positive 
weight  function  w(x)  defined  on  {a,b),  then  an  equidistributed  mesh  is  a  partition 

Iljy  =  =  Xo  <  Atj  <  <  ...  <  x^_j  <  Xfj  =  b)  (2a) 

such  that 

I""'  w(x)dx  =  constant  =  —  f‘‘w(x)dx  ,  j  =  .  (2b) 

Equidistribution  has  been  used  with  functional  approximation  to  determine  the  optimal 
placement  of,  e.g.,  interpolation  points  (cf.,  de  Boor  (ref  33)  and  Soanes  (ref  34)). 
Equidistribution  strategies  have  also  been  used  in  conjunction  with  the  solution  of  two-point 
boundary  value  problems  (cf.,  Pereyra  and  Sewell  (ref  35)).  The  task  of  selecting  a  mesh  to 
minimize  the  discretization  error  is  known  to  be  asymptotically  equivalent  to  equidistributing  the 
local  discretization  error  when  using  the  maximum  norm  (cf.,  de  ^or  (ref  33)). 

Success  with  functional  approximation  and  the  numerical  solution  of  ordinary  differential 
equations  has  led  investigators  to  consider  the  use  of  equidistribution  strategies  for  generating 
moving  meshes  in  conjunction  with  the  numerical  solution  of  transient  partial  differential 
equations  (PDEs)  (cf.,  Da\'is  and  Flaherty  (ref  14)).  The  general  framework  is  to  simply 
reconsider  Eq.  (2)  with  a  time  dependency.  Thus,  the  problem  is  to  determine  a  dynamic  mesh 

n^f)  =  {a  =  Xq  <  x^(f)  <  X2(f)  <  ...  <  Xjy.j(f)  <  x^  =  b)  (3) 

so  that 

w(x,t)dx  -  c(t)  =  f^w(x,t)dx,J  =  ,  t  k  0 ,  (4) 

where  w(x,t)  >  0,  x  e  [a,b\,  r  >  0,  is  usually  chosen  to  be  a  function  of  the  solution  of  the 
underlying  PDE.  For  example,  w  has  been  chosen  to  be  proportional  to  the  solution’s  gradient, 
curvature,  and  local  discretization  error  (cf.,  e.g.,  Adjerid  and  Flaherty  (ref  1),  Bell  and  Shubin 
(ref  13),  Davis  and  Flaherty  (ref  14),  Dorfi  and  Drury  (ref  15),  Dwyer  (ref  16),  Rai  and 
Anderson  (ref  23),  Russell  and  Ren  (ref  24),  Smooke  and  Koszykowski  (ref  26),  Verwer,  Blom, 
Furzeland,  and  Zegeling  (ref  28),  and  White  (ref  29)). 

When  applying  feq.  (4)  in  some  numerical  scheme,  most  investigators  move  a  grid  of  a 
fixed  number  of  elements  to  follow  and  resolve  local  nonuniformities  in  the  solution.  In  order  to 
guarantee  accuracy,  they  must  be  sure  that  N,  the  number  of  elements,  is  large  enough  to 
approximate  the  solution  throughout  the  entire  spatial  domain  and  the  entire  temporal 
integration.  This  is  a  relatively  severe  limitation  since  the  optimal  number  of  elements  is  not 
generally  known  in  advance  and  generally  varies  with  time. 
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Equidistribution  techniques  may  be  applied  to  time-dependent  PDEs  either  by  (1)  solving 
Eq.  (4)  simultaneously  with  the  solution  of  the  PDEs;  (2)  extrapolating  equidistributing  meshes 
at  past  time  levels  to  future  time  levels;  or  (3)  integrating  a  system  of  ordinary  differential 
equations  for,  say,  the  mesh  velocities, 

^[xm  =mj=o,h...,N,  (5) 

at 

that  are  equivalent  to  Eq.  (4).  Many  researchers  have  reported  problems  with  extrapolating 
equidistributing  meshes  or  with  integrating  differential  equations  for  the  mesh  velocities.  For 
example,  if  sufficient  care  is  not  exercised,  mesh  trajectories  can  leave  the  domain  [a,b\,  cross 
each  other,  or  oscillate  wildly  from  time  step  to  time  step.  These  events  can  even  occur  when 
the  solution  of  the  partial  differential  equations  is  changing  very  little. 

Also,  mesh  moving  is  being  used  in  conjunction  with  a  numerical  solution  procedure  for 
the  PDEs.  Since  the  accuracy  of  most  PDE  solution  procedures  depends  on  the  uniformity  of 
the  space -time  grid,  equidistribution  can  deform  the  grid  too  much  and  introduce  new  and  larger 
sources  of  discreti2^tion  errors.  As  a  result,  some  investigators  have  abandoned  mesh  moving  in 
favor  of  local  mesh  refinement  methods  (cf.,  Berger  and  Oliger  (ref  5)  and  Moore  and  Flaherty 
(ref  8)). 

Refinement  of  a  sufficient  level  can  be  used  to  guarantee  that  a  solution  has  a  prescribed 
accuracy;  however,  it  is  more  costly  than  mesh  motion  since  the  solution  must  be  recomputed. 
Additionally,  refinement  is  less  effective  than  mesh  moving  at  reducing  dispersive  errors  in  the 
vicinity  of  wave  fronts. 

Our  choice  has  been  to  combine  local  mesh  refinement  with  mesh  movement  based  on 
equidistribution.  The  benefits  are  twofold.  First,  refinement  will  avoid  any  drastic  deformation 
of  the  grid  caused  by  movement  as  well  as  guarantee  a  prescribed  accuracy.  Second,  mesh 
movement  will  reduce  the  need  for  refinement  and,  thus,  reduce  computational  costs. 

In  the  following  section,  a  method  for  solving  the  static  equidistribution  problem  is 
discussed.  Then  in  the  Mesh  Dynamics:  Linear  Stability  section,  static  equidistribution  is 
extended  to  a  time-dependent  partition,  and  a  linear  perturbation  analysis  is  performed  to 
examine  its  stability.  In  the  Mesh  Dynamics:  Neutral  Stability  section,  dynamic  equidistribution 
is  reformulated  and  its  stability  reexamined  resulting  in  the  construction  of  stable  mesh  moving 
schemes.  Mesh  Movement  Implementation  describes  our  numerical  implementation  of  one  of 
the  stable  mesh  moving  schemes  constructed  in  Mesh  Dynamics:  Neutral  Stability.  Finally,  in  the 
last  section,  a  summary  of  this  report  is  presented. 
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STATIC  SPATIAL  EQUIDISTRIBUTION 


The  equidistribution  problem  Eq.  (2)  can  most  easily  be  solved  by  a  technique  due  to 
de  Boor  (ref  33).  Thus,  we  let 

T{x,t)  =  f'‘w((,t)dC  ,  a  ^  X  ^  b  .  (6) 

J  a 


Then 


c(i)  =  i  nM  (7) 

N 

and  the  equidistributing  mesh  Xj(t),j  =0,1,... ,N,  is  determined  as  the  solution  of  the  nonlinear 
system 

T(Xjit),t)  -  jc(t)  -0,7-  0,1,...^  .  (8) 


The  equidistribution  problem  has  a  nonunique  solution  whenever  w(x,t)  :=  0;  hence,  we 
may  expect  numerical  difficulties  when  w(x,t)  is  small  on  any  subinterval  of  [a,b].  This  problem  is 
usually  handled  by  imposing  a  lower  bound  on  w,  e.g.,  it  is  common  to  replace  w(x,t)  by 
w(x,t)  +  17.  There  are  many  choices  for  the  positive  parameter  rj.  Davis  and  Flaherty  (ref  14) 
suggest  that  77  should  be  related  to  the  discretization  error  of  the  numerical  method  that  is  being 
used  to  solve  the  partial  differential  equations.  Another  popular  choice  (cf.,  e.g.,  Dwyer  (ref  16)) 
is  set  to  77  to  unity,  when  the  interval  [a,b]  and  w  have  been  appropriately  scaled.  Among  other 
things,  both  of  these  choices  ensure  that  the  solution  of  Eq.  (8)  is  a  uniform  mesh  whenever  w  is 
small  everywhere  on  [a,b].  Throughout  the  remainder  of  this  report,  we  assume  that  w(x,t)  >  0 
for  a  <  X  <  6,  t  >  0,  with  w  =  0  only  at  a  finite  number  of  isolated  points.  This  is  sufficient  to 
guarantee  a  unique  solution  of  Eq.  (8). 

If  w(x,t)  is  a  function  of  the  numerical  solution  of  the  associated  partial  differential 
equations,  it  will  generally  only  be  known  discretely.  Suppose  w  is  known  at  the  points  x^,  i  = 

0,1,. ..,M,  then  we  may  approximate  it  by  a  piecewise  polynomial  inx  and  integrate  Eq.  (6)  to  find 
a  piecewise  polynomial  approximation  to  T(x,t).  The  function  c(t)  can  then  be  determined 
approximately  from  Eq.  (7).  An  approximation  7  =  0,1,.. .,N,  to  an  equidistributing  mesh  is 
determined  by  solving  Eq.  (8)  by,  e.g.,  Newton  iteration.  If,  in  particular,  we  approximate  w(x,t) 
by  a  pieeewise  linear  function  of x  with  respect  to  the  mesh  x®,  i  =  0,1,.. .,M,  then  T(x,t)  is  a 
piecewise  quadratic  function  ofx,  and  Eq.  (8)  can  be  solved  for  xj,  j  =  0,1,.. .,N,  directly  by  the 
quadratic  formula  to  give 
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where 


2y  ^ 
P+(p2+2aY)’^  ’ 


j  =  0,1,...^ 


(9a) 


1  0 
Xj  = 


a 


w(XiV)  -  w(x°.i,t) 


(9b) 


P  =  HX°i,t)  y  Y  =J^  -  , 


(9c,d) 


and  i  is  such  that  <  jc  <  T(xfl,t). 

The  number  of  points,  M,  in  the  input  meshx^,  i  =  0,1,.. .,M,  and  N,  in  the  output  mesh 
x^j,  j  =  0,1,...,N,  are  not  necessarily  the  same.  This  could  be  useful  in  situations  where  the 
function  w(x,t)  is  known  very  precisely,  e.g.,  w(x,0)  can  often  be  calculated  exactly  using  the  initial 
data  of  the  associated  partial  differential  equations.  In  this  case,  M  can  be  determined  so  that 
the  integrals  in  Eqs.  (6)  to  (8)  and  the  equidistributing  mesh  can  be  evaluated  to  a  prescribed 
level  of  accuracy.  For  example,  if  the  trapezoidal  rule  is  used  to  approximate  T(x,t),  an 
approximation  of  the  equidistributing  mesh  Xj(t),  j  =  0,1,2,... ,N,  can  be  determined  to  tolerance  e 
by  selecting 

max>v^(x,0 

3^2  >  jb-a)  a<x<b  (10) 

46  min  w(x,t) 

a<x<b 


when  w(x,t)  >  0.  This  estimate  is  obtained  from  standard  error  formulae  for  the  trapezoidal  rule 
(cf.,  Dahlquist,  Bjork,  and  Anderson  (ref  36))  and  elementary  continuity  arguments. 

When  A  =  M,  we  may  think  of  solving  Eqs.  (6)  to  (8)  iteratively.  Thus,  the  mesh  j  = 
0,1,.. .,N,  can  be  used  to  calculate  another  approximations:^,  j  -  0,1,.. .,N,  to  an  equidistributing 
mesh,  etc.  However,  this  iterative  strategy  does  not  necessarily  converge  near  a  local  minimum 
of  w(x,t)  as  illustrated  by  the  following  simple  example. 

Example  1 

Consider  a  three-point  mesh  {M  =  N  =  3),Xj,j  =  0,1,2,  on  -1  <  x  <  1  with  w(x,t)  =  x!^. 
The  endpoints  Xq  =  -1  and  Xj  =  1  are  fixed,  and  the  only  point  that  needs  to  be  determined  by 
iteration  is  Xy  The  exact  value  of  is,  of  course,  zero;  however,  we  start  with  an  initial  guess  x? 
=  e,  use  piecewise  linear  approximations  for  w,  and  see  if  successive  iterates  x],  v  =  1,2,..., 
converge  to  zero.  We  can  show  that: 
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1.  If  e  =  2  :=  2  -  3'^  then =  -z  =  z,v  =  O,!,--  Thus,  the  iterative 

solution  is  a  two-point  limit  cycle. 

2.  If  e  0,  then,  |  <  z,v  =  1,2,.... 

3.  If  <2,  then  K|  < 

4.  If  x\  >  0,  then  x\*'^  <  0. 

Items  (2)  to  (4)  imply  that  jt”  does  not  approach  zero  for  any  nonzero  initial  guess,  but  instead 
approaches  a  limit  cycle,  oscillating  between  z  and  -z  on  alternate  iterations. 


The  illustration  of  above  items  (1)  to  (4)  is  relatively  simple.  The  piecewise  linear 
approximation  to  w(x,t)  =  which  interpolates  to  w  atxj  =  e  is 


-  I  1  ■  ^  -  €),  -1  ^  X  ^  € 

■  \  l  -  (1  -  ;c)(l  +  e),  e  cc  ^  1  • 

(11) 

This  implies  that 

71:€,0  =  ^(1  +  £)(1  +  £2), 

(12a) 

7K1,?)  =  1  +  £2, 

(12b) 

and 

c(I)  =  j(l  +  e"). 

(12c) 

Now  consider  e  >  0  (the  case  e  <  0  follows  from  symmetry),  then 

7X0,0  =  |(1  +  6). 

(13) 

Since  0  <  e  <  1,  c(t)  <  T(0,t),  implying  <  0.  This  proves  item  (4). 

The  solutions}  is  such  that  T(x^i,t)  =  c(t),  i.e., 

(1  +  xl)  -  ^(1  +  xl)Hl  -  e)  =  ^(1  +  e^).  (14) 


This  implies  that 
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1  ^  6  -  Vl  -  (l-£)(l+£^) 
1  -  € 


^1  = 


(15) 


Now  one  might  ask  if  there  exists  an  e  such  that  x\  =  -e  7  The  answer  is  found  by 
substituting  -e  for  in  Eq.  (14).  This  substitution  yields 


(1  -  £)  -  i(l  -  e)’  =  i(l  +  e^). 
2  2 


(16a) 


or 


€  -  4e^  +  =  0. 


(16b) 


Equation  (16)  factors  into  the  following  equation: 

e[e  -  (2-v^)][€  -  (2+^3)]  =  0. 

Since  z  :=  2  -  3''""  is  the  only  factor  in  the  open  interval  (0,1),  this  result  proves  item  (1). 
Equation  (16c)  implies  that  if  0  <  e  <  z,  then 

€[€  -  (2-^/3)][€  -  (2+V^)]  >  0. 


(17) 


Equation  (17)  implies  that  T(-€,t)  >  c(t).  Therefore,  x}  <  -e  whenever  0  <  e  <  z.  This  proves 
item  (3). 

All  that  remains  is  to  prove  item  (2).  Rewrite  e  as  z  +  e^,  where  is  such  that  -1  <  z  + 
<  1.  Then  item  (2)  is  true  if  T(-z,t)  <  c(t).  This  implies  that  item  (2)  is  true  if 


(1  -  Z)  -  |(1  -  Z)'‘(l  -  (Z-€^))  <  |(1  +  (Z+€/).  ' 


(18) 


Equation  (18)  reduces  to 


|(1  -2z+Z^)e,  <  |(2z€^  +  eJ) 


(19a) 


or 


0  <  -|(1 -4z+z2)€^  + 

Since  z  is  a  root  of  the  quadratie  term  of  Eq.  (19b),  this  implies  that  item  (2)  is  true  if 


(19b) 
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(20) 


0  < 


This  proves  item  (2). 


MESH  DYNAMICS:  LINEAR  STABILITY 


The  discussion  of  the  previous  section  involved  the  computation  of  an  equidistributing 
mesh  at  one  time.  As  noted  in  the  Introduction,  equidistributing  meshes  at  subsequent  time 
levels  can  be  obtained  by  a  variety  of  techniques.  In  this  section,  we  use  linear  stability  analyses 
to  explain  why  many  researchers  have  been  experiencing  difficulties  when  using  either 
extrapolation  or  ordinary  differential  equations  for  mesh  velocities  to  obtain  equidistributing 
meshes.  There  are  no  stability  problems  associated  with  solving  Eq.  (8)  directly;  however  an 
objection  to  this  approach  is  the  complexity  associated  with  solving  differential-algebraic  systems 
for  the  partial  differential  equation  solution  and  the  mesh  motion. 


We  begin  our  analysis  by  considering  a  system  obtained  by  differentiating  Eq.  (8)  with 
respect  to  time.  Upon  use  of  Eq.  (6),  this  gives 


w(Xj,t)Xj  =  -J‘\p,f)dx  -  jc(t) 


,  j  =  u,...^-l 


(21) 


Suppose  Xj(t),  j  =  0,1,...,N,  is  an  equidistributing  mesh  that  exactly  satisfies  Eqs.  (8)  and 
(21)  for  t  >  0.  Let  us  introduce  a  small  perturbation  Sxj(O),  j  =  0,1,. ..,N,  at  r  =  0  and  study  its 
effect  on  Eq.  (21).  If  no  additional  errors  are  introduced,  the  perturbed  system  satisfies 

w{xft)  +bxft),t)  ix.+hxp  =  - ^J^'*^\ix,t)dx  -jc(t)  ,  (22) 

subject  to  the  constraints  Sxjt)  =  Sxjt)  =  0.  We  further  assume  that  <<  b-  a,  j  = 

1, 2,.. .,N-1,  and  linearize  Eq.  (22)  to  get 

■^Mx/t),t)bxj(t)]  =0,7  =  .  (23) 

Integrating,  we  find 

bXj(f) 


W(Xj(t),t) 


bx/0)  ,j  =  l,2,...,Af-l 


(24) 


Therefore,  the  differential  system  Eq.  (21)  is  stable  to  linear  perturbations  if 
is  less  than  unity.  Unfortunately,  most  choices  of  w(x,t)  are  likely  to  be  decaying  functions  of 
time  for  dissipative  parabolic  systems,  and  this  will  almost  certainly  yield  a  value  of  L(t)  that  is 
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(25) 


max  ^■*y((^)>0) 
Mxft),t) 


larger  than  unity.  Local  instabilities  can  also  occur  whenever  the  mesh  is  moved  so  that  L(t) 
exceeds  unity  for  some  specific  times.  These  instabilities  may  grow  or  decay  as  time  progresses 
depending  on  the  value  of  L. 

The  following  two  examples  illustrate  some  of  the  instabilities  that  can  occur. 

Example  2 


Consider  the  heat  conduction  problem 


(26a) 

u(x,0)  =  sinnx  ,  u(0,t)  =  u(l,t)  =  0  . 

(26b, c) 

The  exact  solution  of  this  problem  is 

u(x,t)  =  e  ‘sinnx  . 

(26d) 

We  take 

H<X,f)  «  |M^(x,0i  . 

(26e) 

Since  u(x,t)  and  w(x,t)  are  separable  in  space  and  time,  the  correct  strategy  is  to  generate  an 
equidistributed  mesh  at  time  t  =  0  and  use  it  for  all  time.  However,  L(t)  «  exp(-t),  and  thus  we 
expect  the  solution  of  Eq.  (21)  to  be  unstable.  In  Figure  1  (at  the  end  of  this  section),  we 
display  the  meshes  produced  by  both  Eqs.  (8)  (dashed  curves)  and  (21)  (solid  curves),  and  the 
unstable  behavior  of  Eq.  (21)  is  clearly  visible.  Simpson’s  rule  with  M  =  200  was  used  to 
evaluate  all  integrals,  a  mesh  of  TV  =  10  elements  was  equidistributed,  and  an  initial  perturbation 
Sxj(O)  =  0.025,  j  =  i,2,...,TV-i,  was  introduced.  Equations  (21)  (and  all  differential  equations 
appearing  in  Example  3)  were  solved  using  a  version  of  the  differential-algebraic  equation  solver 
DASSL  (cf.,  Petzold  (ref  37)). 

Example  3 

We  consider  a  problem  for  a  partial  differential  equation  that  has  the  exact  solution 

u(x,t)  =  tanhrj[(x-l)+r2f]  ,0sx^l,f^0.  (27) 

The  function  (27)  is  a  wave  that  travels  in  the  negative  x  direction  when  and  rj  are  positive  (cf.. 
Figure  2  at  the  end  of  this  section).  The  values  of  and  ^2  determine  the  steepness  of  the  wave 


9 


and  its  propagation  speed.  Such  solutions  could  arise  from  several  types  of  partial  differential 
equations,  e.g.,  Davis  and  Flaherty  (ref  14)  studied  a  heat  conduction  problem  of  the  form 

+  Axit)  =  \  ^  X  ^  \  ,t>  0  ,  (28) 

<r 

where  the  initial  conditions,  Dirichlet  boundary  conditions,  constant  diffusion  ct and  source  / 
were  chosen  so  that  the  exact  solution  was  given  by  Eq.  (27). 

The  meshes  produced  by  both  Eqs.  (8)  (dashed  curves)  and  (21)  (solid  curves)  for  =  5, 
r2  =  1  and 

Mx,t)  oc  \uj^x,t)\  +  q  ,  (29) 

where  tj  =  0.1,  are  shown  in  Figure  3  (at  the  end  of  this  section).  The  solution  of  Eq.  (21)  is 
marginally  stable  for  small  times,  but  w  decreases  as  the  wave  front  progresses  towards  x  =  0, 
and  the  instability  is  apparent.  In  fact,  some  mesh  trajectories  have  left  the  domain  [0,1]. 
Simpson’s  rule  with  M  =  200  was  used  to  evaluate  all  integrals,  a  mesh  of  =  10  elements  was 
equidistributed,  and  an  initial  perturbation  8x^(0)  =  0.025,  j  =  1,2,...,N-1,  was  introduced. 
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1.8 


ctories  for  Example  3. 
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MESH  DYNAMICS:  NEUTRAL  STABILITY 


In  the  previous  section,  it  was  shown  how  in  certain  instances  mesh  moving  could  become 
unstable.  Such  instabilities  help  to  explain  the  oscillations  that  some  researchers  have  observed 
with  mesh  movement.  However,  other  researchers  have  seen  no  such  instabilities  develop  in 
their  mesh  moving  schemes,  even  when  solving  dissipative  problems.  The  results  of  the  previous 
section  do  not  give  sufficient  guidance  to  support  the  construction  of  stable  mesh  moving 
schemes.  In  order  to  explore  this  further  and  to  more  fully  understand  the  dynamics  of  mesh 
movement,  define 


TTxfl  /VC.<)rfC 

=  if - 

nb,t) 

J  a 

(cf.,  Eq.  (6))  and  rewrite  Eq.  (8)  as 

^(x(t),t)  ~  J-  =  0  J  =  0,1.... JV. 
N 

Now  consider  the  differential  system, 

^  0 ,  i  ^ 

at  ' 


(30) 


(31) 


(32) 


obtained  by  differentiating  Eq.  (31)  with  respect  to  time. 


The  next  step  is  to  perform  the  same  stability  analysis  on  Eq.  (32)  as  was  done  on  Eq. 
(21)  in  the  previous  section.  This  perturbation  analysis  can  be  simplified  by  noting  that  Eq.  (31) 
is  identical  to  Eq.  (8)  with  the  weight  function  w{x,t)  replaced  by  the  normalized  weight  function 


W(x,t)  = 


Hx,f) 

J  a 


(33) 


Recall  that  w(x,t)  is  zero  at  only  a  finite  number  of  isolated  points;  hence,  W(x,t)  is  well  defined. 
Therefore,  the  result  of  the  stability  analysis  for  Eq.  (32)  can  be  found  by  replacing  w  by  IE  in 
Eq.  (25).  The  differential  system  Eq.  (32)  is  stable  to  linear  perturbation,  then,  if 


max  ^(^^(0)>0) 


(34) 


is  less  than  unity. 
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At  first  glance,  it  would  appear  that  the  stability  determination  has  simply  shifted  from 
w(x,t)  to  W(x,t);  thus,  instability  occurs  when  L(t)  >  1  or  when  W{x-{t),t)  is  a  decreasing  function 
of  time.  However,  if  Eq.  (34)  is  rewritten  with  (33)  used  to  indicate  the  dependence  of 
w(x,t),  the  result  is 


max  vK;t/o),o)  /Vw; 

0<j<N 

J  a 


(35) 


The  conditions  expressed  by  Eq.  (35)  are  potentially  less  severe  than  those  of  Eq.  (25)  since 


m  ^  m 


(36) 


whenever  w(x,t)  is  a  decreasing  function  in  time.  Some  improvement  may  have  occurred, 
although  there  is  still  no  guarantee  of  stability.  Significant  improvement  does  occur  when  w(x,t) 
is  a  separable  function  of  space  and  time  (cf..  Example  2).  In  this  case,  L(t)  =  1;  thus,  mesh 
motion  is  neutrally  stable  whether  w(x,t)  h  decreasing  or  not.  Thus,  for  any  separable  function 
w(x,t),  small  perturbations  of  Eq.  (32)  will  not  grow  or  decay  as  time  increases. 

Example  4 

Reconsider  Example  2  where  w(x,t)  was  given  by  Eq.  (26).  Recall  that  Eq.  (21) 
generated  unstable  mesh  motion  because  L(t)  «  exp(-t).  Since  L(l)  =  1  for  this  function  w(x,t), 
the  solution  of  Eq.  (32)  should  show  no  evidence  of  these  instabilities.  In  Figure  4  (at  the  end  of 
this  section)  meshes  produced  by  both  Eqs.  (31)  (dashed  curves)  and  (32)  (solid  curves)  are 
displayed.  As  predicted,  mesh  trajectories  generated  by  Eq.  (32)  remain  neutrally  stable  for  all 
time.  As  in  Example  2,  Simpson’s  rule  with  M  =  200  was  used  to  evaluate  all  integrals,  a  mesh 
of  N  =  10  elements  was  equidistributed,  and  an  initial  perturbation  of  8xj(0)  =  0.025,  j  =  1,...,10, 
was  introduced.  Equation  (32)  was  solved  using  a  version  of  the  differential-algebraic  solver, 
DASSL  (cf.,  Petzold  (ref  37)). 

Equidistribution  through  Eq.  (32)  does  not  guarantee  that  initial  small  perturbations  do 
not  grow  unless  w(x,t)  is  separable.  However,  if  Eq.  (32)  is  reexamined,  with  emphasis  placed  on 
the  stability  of  the  functions  ^{xj(t),t),  j  =  the  results  are  different.  In  order  to  see  this, 

define 

=  <I>(:c/f),0  ,j  =  (37) 

and  rewrite  Eq.  (32)  as 


=0,7  =  .  (38) 

dt  ^ 

The  stability  results  are  obvious  in  terms  of  these  new  dependent  variables.  Equation  (38)  is 
neutrally  stable  and  initial  perturbations  of  (through  x^,j  =  will  neither  grow  nor 

decay.  This  does  not  imply  that  the  mesh  trajectories  are  neutrally  stable.  However,  a  type  of 
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(39a) 


stability  is  inherited.  If  Sx/O),  j  =  1,2,...,N-1,  is  such  that 

<  x/0)  +  6x/0)  <  x.,j(0)  ,  j  =  IX-.  Jf-l  , 

then  solutions  of  Eq.  (33)  will  ensure  that 

<  xft)  +  bxft)  <  ,  j  =  IX-.  Jf-l  ,  (39b) 

for  all  t  >  0. 

Example  5 

Reconsider  Example  3  where  w(x,t)  was  given  by  Eq.  (29).  Recall  that  Eq.  (21) 
generated  unstable  mesh  motion  because  w  decreased  as  the  wave  front  progressed  to  x  =  0. 
Solutions  of  Eq.  (38)  should  demonstrate  no  such  instabilities.  In  Figure  5  (at  the  end  of  this 
seetion)  meshes  produced  by  both  Eqs.  (31)  (dashed  curves)  and  (38)  (solid  curves)  are 
displayed.  As  predicted,  mesh  trajectories  generated  by  Eq.  (38)  remain  bounded  for  all  time 
and  do  not  intersect.  As  in  Example  4,  Simpson’s  rule  with  M  =  200  was  used  to  evaluate  all 
integrals,  a  mesh  of  N  =  10  elements  was  equidistributed,  and  an  initial  perturbation  of 
&Xj(0)  =  0.025,  j  =  1,2,. ..,10,  was  introduced.  Equation  (38)  was  solved  using  a  version  of  the 
differential-algebraic  solver,  BASSE  (cf.,  Petzold  (ref  37)). 

Furthermore,  if  Eq.  (21)  is  rewritten  in  terms  of  4>, ,  j  =  1,2,...,N-1,  the  result  is 

‘  i-MD  ,J  ‘  1X-.N-1  ,  (40a) 

at  ^  ^  N 

where 

[\iU)dC 

A(t)  =  ^ -  .  (40b) 

J  a 


Equation  (40)  will  be  asymptotically  stable  if  A(t)  >  1  and  neutrally  stable  when  A(t)  =  1.  Initial 
perturbations  will  grow  when  A(t)  <  1.  The  condition,  A(t)  <  1,  will  occur  when  w,(x,t)  <  0,  i.e., 
when  w(x,t)  is  a  decreasing  function  of  time.  This  result  is  stronger  than  that  of  the  previous 
section,  since  arbitrary  perturbations  are  permitted. 

Equation  (40)  also  demonstrates  how  to  form  asymptotically  stable  mesh  moving 
equations,  if  so  desired.  Simply  substitute  a  positive  constant,  >  0,  for  A(t)  in  Eq.  (40a)  to 
yield 


at  ^  ^  N 


(41) 


Now,  not  only  are  solutions  to  Eq.  (41)  stable,  but  can  be  chosen  to  adjust  the  decay  rate  of  a 
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perturbation  to  its  equidistributed  position.  This  idea  is  similar  to  one  proposed  by  Hyman  (cf., 
Coyle,  Flaherty,  and  Ludwig  (ref  38)).  Notice  that  unlike  the  neutrally  stable  case  (cf.,  Eq.  (38)), 
an  initial  perturbation  of  ,  j  =  l,2,...,N-ly  decaying  to  zero  does  imply  that  the  associated 
initial  mesh  perturbation  of Xj,j  =  i, 2, ...,N-i, decays  to  zero. 

Example  6 

Reconsider  Example  4  where  w(x,t)  was  given  by  Eq.  (26).  Recall  that  Eq.  (32) 
generated  neutrally  stable  mesh  motion  because  L(t)  =  1.  Solutions  of  Eq.  (41)  should 
asymptotically  approach  the  equidistributed  positions.  In  Figure  6  (at  the  end  of  this  section), 
meshes  produced  by  both  Eqs.  (31)  (dashed  curves)  and  (41)  (solid  curves)  are  displayed.  As 
predicted,  mesh  trajectories  generated  by  Eq.  (41)  approach  those  generated  by  Eq.  (31)  as  time 
progresses.  As  in  feample  4,  Simpson’s  rule  with  M  =  200  was  used  to  evaluate  all  integrals,  a 
mesh  ofN=  10  elements  was  equidistributed,  and  an  initial  perturbation  of  Sxj(O)  =  0.025,  j  = 
1,2,...,N-1,  was  introduced.  A  value  of  =  1  was  used,  and  Eq.  (41)  was  solved  using  a  version 
of  the  differential-algebraic  solver,  DASSL  (cf.,  Petzold  (ref  37)). 

Example  7 

Reconsider  Example  5  where  w(x,t)  was  given  by  Eq.  (29).  Recall  that  Eq.  (38) 
generated  neutrally  stable  mesh  motion.  Solutions  of  Eq.  (41)  should  asymptotically  approach 
the  equidistributed  position.  In  Figure  7  (at  the  end  of  this  section),  meshes  produced  by  both 
Eqs.  (31)  (dashed  curves)  and  (41)  (solid  curves)  are  displayed.  As  predicted,  mesh  trajectories 
generated  by  Eq.  (41)  approach  those  generated  by  Eq.  (31)  as  time  progresses.  As  in  Example 
5,  Simpson’s  rule  with  M  =  200  was  used  to  evaluate  all  integrals,  a  mesh  ofN=  10  elements 
was  equidistributed,  and  an  initial  perturbation  of  Sx/O)  =  0.025,  j  =  i,2,...,N-i,was  introduced. 

A  value  of  A^  =  1  was  used,  and  Eq.  (41)  was  solved  using  a  version  of  the  differential  algebraic 
solver  DASSL  (cf.,  Petzold  (ref  37)). 

Another  positive  feature  of  these  schemes  is  that  extending  them  to  higher  dimensions  is 
easily  accomplished  with  all  of  their  stability  properties  left  intact.  For  simplicity,  consider  a 
rectangular  two-dimensional  domain 

fl  =  {(A:,y)|a  ^x^b,c^y^d} 

and  the  extension  of  Eq.  (38)  to  SI.  To  that  end,  define 

f^w(C,r],t)dr)dC 

T(W)  =  -  (43) 

/7  VC,Ti,0rfTirfC 

J  c  ^  c 

where  w(x,y,t)  is  a  positive  weight  function  on  SI  such  that  ^^(x,y,t)  is  a  well-defined  function  with 
the  necessary  continuity  properties.  Furthermore,  let 
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(44a) 


IL\t)  =  {(xp),y^(t))}  ,j  =  0,1,...^^  ,  k  = 

denote  a  partition  of  ft  into  subrectangles  at  any  time  t  such  that 

a  =  x^it)  <  x^(t)  <  ...  <  x^(t)  =  b  , 


(44b) 


c  =  yoit)  <  )’i(r)  <  yjjjt)  =  d  .  (44c) 

In  order  for  n^(t)  to  move  in  a  neutrally  stable  manner  for  any  initial  partition  n^(0), 
demand  that  it  obey  the  following  equations  for  r  >  0; 


dt  ^ 

(45a) 

jmb,y^(tm  =0,7  =  l,2,...,Afy-l  . 

(45b) 

Note  that  the  stability  of  the  two-dimensional  movement  is  guaranteed  by  the  one-dimensional 
theory  since  Eq.  (45a)  is  equivalent  to  one-dimensional  movement  in  the  j;-direction,  while  Eq. 
(45b)  is  equivalent  to  one-dimensional  movement  in  the  y-direction. 
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MESH  MOVEMENT  IMPLEMENTATION 


Recall  that  the  purpose  of  investigating  mesh  moving  schemes  was  to  develop  an  effective 
tool  to  aid  in  the  process  of  solving  PDEs  numerically.  In  the  Mesh  Dynamics:  Linear  Stability 
and  Mesh  Dynamics:  Neutral  Stability  sections,  the  stability  of  mesh  movement  was  discussed  in 
isolation  from  its  use  with  a  PDE  solver.  The  reason  was  that  if  mesh  movement  is  not  stable, 
then  there  would  be  no  hope  of  it  being  a  useful  numerical  tool.  Now  that  the  stability  question 
has  been  answered,  the  question  of  implementation  in  a  PDE  solver  needs  to  be  addressed. 

The  first  consideration  is  the  choice  of  a  mesh  moving  scheme.  Since  an  unstable  scheme 
is  not  an  option,  the  choice  would  seem  to  be  between  a  neutrally  stable  scheme  (cf.,  Eq.  (38)) 
or  an  asymptotically  stable  one  (cf.,  Eq.  (41)). 

Usually  in  the  field  of  numerical  analysis,  asymptotically  stable  schemes  are  preferred  (cf., 
e.g.,  Lapidus  and  Pinder  (ref  39)  and  Richtmyer  and  Morton  (ref  40)).  This  is  because 
computational  errors  are  unavoidable,  and  asymptotically  stable  methods  are  best  at  minimizing 
the  consequences  of  those  errors.  However,  in  this  instance,  the  choice  of  a  mesh  moving 
scheme  should  be  more  related  to  the  nature  of  the  associated  PDE  solver. 

If  the  solver  is  topologically  invariant,  i.e.,  mesh  refinement  is  not  present  (cf.,  Davis  and 
Flaherty  (ref  14)  and  Miller  and  Miller  (refs  20,21)),  then  an  asymptotically  stable  mesh  moving 
scheme  is  indeed  the  method  of  choice.  Since  the  dimension  of  the  mesh  is  fixed,  the  best  one 
can  do  is  position  them  optimally.  This  is  precisely  the  situation  where  an  asymptotically  stable 
scheme  works  best. 


However,  if,  as  is  proposed  here,  refinement  will  be  performed  in  conjunction  with  mesh 
movement,  then  a  neutrally  stable  scheme  is  preferred.  This  is  because  a  neutrally  stable  scheme 
is  independent  of  the  number  of  mesh  points  (cf.,  Eq.  (38)).  Also,  there  is  no  unknown 
parameter  A  to  choose. 


Suppose,  for  example,  that  refinement  demands  a  new  point,  to  be  added  at  time  f 
between  existing  points  xj  and  To  determine  how should  move  for  t  >  f,  simply 
solve  the  following  initial  value  problem: 


(46a) 
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If,  on  the  other  hand,  refinement  requires  deletion  of  a  point,  then  mesh  motion  continues  with 
the  remaining  points  moving  according  to  Eq.  (38). 
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In  this  way,  theoretically,  when  assuming  exact  integration,  the  addition  or  deletion  of 
points  via  refinement  has  no  effect  on  the  movement  of  any  previously  existing  points  that 
remain  after  refinement.  The  movement  of  new  points  will  begin  at  the  instant  of  their  creation 
and  will  be  unaffected  by  anything  the  previously  existing  points  have  been  or  will  be  doing.  As 
a  consequence,  we  chose  to  implement  the  neutrally  stable  scheme  Eq.  (38)  over  the 
asymptotically  stable  one. 

The  next  consideration  is  the  choice  of  the  equidistribution  function  w(x,t)  (cf.,  Eq.  (30)). 
This  choice  should  be  related  to  the  reasons  for  incorporating  mesh  movement  into  the  PDE 
solver  since  w(x,t)  will  determine  the  qualitative  nature  of  the  movement  (cf.,  Eq.  (38)  or  (41)). 

In  general,  the  aim  of  any  adaptive  strategy  is  to  increase  the  reliability  and  reduce  the 
temporal  and  spatial  complexity  of  the  solution  procedure.  With  h-refinement,  this  takes  the 
form  of  solving  a  problem  to  some  specified  accuracy  using  a  coarser  discretization  level  than 
would  be  necessary  without  adaptivity. 

There  seems  to  be  no  general  consensus,  however,  as  to  what  w(x,t)  should  be  in  order  to 
simplify  the  solution  process.  Some  researchers,  for  example,  choose  w(x,t)  to  be  proportional  to 
some  physical  quantity  of  interest  and/or  proportional  to  some  derivative  of  the  solution  (cf..  Bell 
and  Shubin  (ref  13),  Davis  and  Flaherty  (ref  14),  Dorfi  and  Drury  (ref  15),  Dwyer  (ref  16),  Rai 
and  Anderson  (ref  23),  Russell  and  Ren  (ref  24),  Smooke  and  Kos^kowski  (ref  26),  Verwer, 
Blom,  Furzeland,  and  Zegeling  (ref  28),  and  White  (ref  29)).  Still  others  attempt  to  relate  w(x,t) 
to  the  residual  error  (cf..  Miller  and  Miller  (refs  20,21))  or  the  local  discretization  error  (cf., 
Adjerid  and  Flaherty  (ref  1)).  The  underlying  assumption  with  any  choice  is  that  reasonable 
approximations  to  the  chosen  monitors  exist. 

When  error  estimates  are  available  (cf.,  Adjerid  and  Flaherty  (ref  1)),  it  would  seem 
appropriate  to  let  these  estimates  control  the  movement.  However,  more  often  than  not,  these 
estimates  lead  to  erratic  mesh  motion  due  to  their  dependency  on  mesh  location.  Some 
researchers  reduce  this  difficulty  by  coupling  the  error  estimation  and  mesh  moving  process  with 
the  PDE  solution  (cf.,  Adjerid  and  Flaherty  (ref  1)). 

Experience  has  shown  that  the  estimates  developed  for  use  in  AFEM  (cf.,  Coyle  and 
Flaherty  (ref  30))  do  not  produce  reliable  mesh  motion,  when  left  uncoupled  from  the  solution 
process.  However,  our  desire  is  to  develop  a  mesh  moving  strategy  that  is  separate  from  the 
PDE  solver,  not  because  coupling  is  ineffective  (cf.,  Adjerid  and  Flaherty  (ref  1)),  but  rather 
because  of  the  greater  flexibility  provided  by  an  uncoupled  approach  to  mesh  moving.  In  this 
way,  movement  can  be  simple  to  implement,  easily  invoked  or  ignored  during  the  solution 
process,  transparently  combined  with  a  variety  of  PDE  solvers  which  have  no  mesh  movement 
capabilities,  and  extensible  to  higher  dimensional  problems. 

When  error  estimates  are  unavailable  or  inappropriate,  picking  w(x,t)  to  be  proportional 
to  some  combination  of  the  approximate  solution  and/or  some  of  its  derivatives  seems  plausible. 
This  is  because  the  discretization  error  of  any  PDE  solver  can  usually  be  expressed  in  terms  of 
some  solution  derivatives  (cf.,  e.g.,  Lapidus  and  Pinder  (ref  39),  Richtmyer  and  Morton  (ref  40), 
and  Strang  and  Fix  (ref  41)). 
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At  first,  we  proceeded  along  lines  similar  to  Davis  and  Raherty  (ref  14),  and  chose  w  to 
be  proportional  to  the  second  spatial  derivative  of  the  approximate  solution.  The  spatial  error  of 
a  piecewise  linear  approximation  is  proportional  to  the  corresponding  exact  spatial  derivative  so 
equidistributing  such  a  w  should  reduce  the  spatial  error  component  of  the  method.  However, 
when  left  uncoupled  from  the  solution,  w  must  be  determined  at  an  advanced  time  level  from 
past  data  in  order  to  propagate  the  mesh  to  that  advanced  time  level.  Experiments  showed  that 
although  this  strategy  was  successful  for  the  first  few  time  steps,  numerical  errors  accumulated  to 
the  point  where  the  extrapolation  of  this  w  would  produce  hi^ly  inaccurate  approximations  to 
the  second  spatial  derivative  and,  hence,  unreliable  mesh  motion. 

Experiments  using  the  exact  second  spatial  derivative  for  problems  where  the  exact 
solution  was  known  proved  revealing.  Although  the  movement  did  reduce  the  spatial  error 
component,  the  temporal  error  component  of  the  method  increased  to  the  point  where  the  total 
error  actually  increased  as  well.  This  was  unacceptable.  As  a  result,  we  began  to  investigate 
choices  for  w  that  addressed  temporal  as  opposed  to  spatial  concerns  (cf.,  e.g.,  Ewing,  Russell, 
and  Wheeler  (ref  17),  Hyman  (ref  18),  Kansa,  Morgan,  and  Morris  (ref  19),  and  Petzold  (ref 
22)).  This  also  had  the  appeal  of  simplifying  the  extension  to  higher  dimensions  because  only  the 
spatial  dimensions  increase  and  the  nature  of  the  temporal  dimension  remains  the  same. 

Choosing  w  to  be  proportional  to  first  and/or  second  temporal  derivatives  of  the 
approximate  solution  proved  unsatisfactory.  Again,  experiments  indicated  that  extrapolation 
would  eventually  produce  inaccurate  approximations  to  these  temporal  derivatives. 

Experimentation  with  exact  temporal  derivatives  was  successful,  however.  The  temporal 
component  of  the  error  decreased  and  the  spatial  component  did  not  increase  significantly,  so 
the  total  error  also  decreased  as  a  result  of  the  movement  (cf.,  Coyle  (ref  42)). 

Only  when  extrapolating  the  approximate  solution  to  advanced  time  levels  were  numerical 
errors  small  enough  so  as  not  to  degrade  the  approximation  significantly.  This  is  because  the 
approximate  solution  is  known  to  a  higher  order  of  accuracy  than  its  derivatives. 

At  first,  choosing  w  to  be  proportional  to  the  approximate  solution  did  nof  seem  so 
appealing  since  discretization  errors  are  more  closely  related  to  derivatives  of  solutions. 

However,  other  researchers  have  had  success  when  moving  the  mesh  so  that  the  variation  in  time 
of  the  solution  is  reduced  in  the  transformed  coordinates  (cf.,  e.g.,  Ewing,  Russell,  and  Wheeler 
(ref  17),  Hyman  (ref  18),  Kansa,  Morgan,  and  Morris  (ref  19),  and  Petzold  (ref  22)).  This  seems 
plausible  since  a  slowly  varying  solution  implies  small  temporal  derivatives. 

Upon  reexamining  Eq.  (38),  one  sees  that  in  the  transformed  coordinates,  the  quantities 
<J>,(t),;  =  are  demanded  to  remain  constant  not  just  slowly  varying.  Such  a  scheme, 

with  w  chosen  to  be  proportional  to  the  approximate  solution,  should  also  keep  temporal 
derivatives  small,  and,  hence,  temporal  errors,  as  experimentation  has  testified  (cf.,  Coyle  and 
Flaherty  (ref  32)).  In  addition,  when  considering  the  extension  to  higher  dimensions,  no  extra 
derivatives  enter  the  computation. 

Our  decision  has  been  to  let  w(x,t)  be  proportional  to  the  L2  norm  of  the  approximate 
solution.  Since  Eq.  (38)  simply  states  to  keep  <I>,(t),y  =  i, 2, ...,/V-7,  constant  from  time  step  to 
time  step,  the  following  procedure  is  performed.  First,  linearly  extrapolate  the  mesh  to  time 
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level  from  mesh  positions  at  time  levels  f  and  If  mesh  trajectories  cross  or  leave  the 
domain,  then  delete  these  offending  points.  This  first  mesh  is  not  intended  to  be  the 
computational  mesh  but  rather  a  mesh  on  which  to  determine  the  extrajX)lated  values  of 
w(x,f'^^).  As  such,  nothing  sophisticated  is  performed  to  avoid  possible  mesh  instabilities. 

Next,  linearly  extrapolate  the  approximate  solution  onto  the  mesh  at  time  level  from 
data  at  time  levels  r  and  Then  determine  w(x,f^‘)  as  the  L2  norm  of  this  extrapolated 
solution. 

Finally,  the  new  mesh  points  at  time  level  f*’  are  determined  such  that 

*/(”•')  -*/»'),;■-  1, 2, .... W-l  .  (47) 

Eq.  (47)  is  solved  by  the  method  detailed  in  the  Static  Spatial  Equidistribution  section. 


SUMMARY 

An  extension  of  one-dimensional,  static  mesh  equidistribution  to  a  time-dependent 
domain  was  proposed.  The  purpose  of  this  extension  was  to  provide  aid  in  the  process  of  solving 
partial  differential  equations  numerically. 

The  stability  of  this  extension  was  examined.  As  a  result,  we  explained  why  many 
intuitively  obvious  schemes  for  calculating  equidistributing  meshes  for  time-dependent  partial 
differential  equations  are  unstable.  In  particular,  we  derived  a  stability  condition  on  the 
equidistribution  density  w(x,t)  that  can  easily  be  verified  in  practice.  In  addition,  we  constructed 
both  neutrally  stable  and  asymptotically  stable  mesh  moving  techniques  as  demonstrated  by  the 
various  examples  of  the  Mesh  Dynamics;  Neutral  Stability  section,  as  well  as  projX)sed  an 
extension  to  higher  dimensional  domains.  The  neutrally  stable  scheme  was  incorporated  into  a 
numerical  method  for  solving  PDEs  that  already  possessed  error  estimation  and  mesh  refinement 
capabilities. 

The  results  of  this  report  indicate  that  mesh  movement  is  still  a  credible  adaptive  option. 
The  problems  researchers  have  experienced  with  mesh  moving  are  not  because  of  any  intrinsic 
failure  of  the  approach,  but  rather,  are  a  result  of  a  lack  of  understanding  of  the  stability 
properties  of  the  process.  This  led  to  the  stability  analyses  of  Mesh  Dynamics:  Linear  Stability 
and  Mesh  Dynamics;  Neutral  Stability  sections  as  well  as  to  the  construction  of  stable  moving 
schemes. 

Mesh  movement  cannot  guarantee  that  a  computed  solution  is  accurate  to  some 
prescribed  tolerance,  however,  because  the  discretization  level  is  fixed.  In  order  to  overcome  this 
obstacle,  a  local  mesh  refinement  algorithm  with  error  estimation  capabilities  has  also  been 
incorporated  into  AFEM.  Local  mesh  refinement  schemes  can  be  costly,  due  to  the  necessity  of 
recomputing  the  solution;  however,  proper  mesh  motion  should  permit  as  few  levels  of 
refinement  as  possible  (cf.,  Coyle  and  Flaherty  (ref  32)). 
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